#define USE_PUPY_WEIGHTS true
#define MIN_GOOD_SCORE 10.0

#include "nw_aligner_simple.cpp"
#include "nw_aligner_interpolation.cpp"

int main( int argc, char *argv[] ) {
    string input_line, genome, genome_id, gene, gene_id, genome_subseq, prev_gene_id;
    vector<string> splited_str;
    double score, raw_score, prev_score;
    int start_pos, length;

    const char *interpolationFileName = NULL;    /* параметр -i */
    static const char *optString = "i:";
    bool use_external_interpolation = false;
    InterpolationData interpolation_data;
    InterCoeffsMC ikoeffs;
    prev_gene_id = "";
    prev_score = 0.0;

    // read options
    int opt = 0;
    opt = getopt( argc, argv, optString );
    while( opt != -1 ) {
        switch( opt ) {
            case 'i':
                interpolationFileName = optarg;
                use_external_interpolation = true;
                break;
            default:
                /* сюда на самом деле попасть невозможно. */
                break;
        }
        opt = getopt( argc, argv, optString );
    }

    if (use_external_interpolation) {
        if (0 != init_interpolation_data(interpolation_data, interpolationFileName)) {
            std::cerr << "Unable to init from interpolation file \"" << interpolationFileName << "\""<< std::endl;
            return 1;
        }
        else {
            //std::cerr << "Use external interpolation from file \"" << interpolationFileName << "\""<< std::endl;
        }
    }
    else {
        std::cerr << "Use internal Monte-Carlo" << std::endl;
    }


    getline(cin, input_line);
    if (input_line.size() == 0) {
        return 0;
    }
    
    splited_str = split(input_line, '\t');
    if (splited_str.size() != 2) {
        cerr << "Bad genome line!" << endl;
        return 1;
    }
    genome_id = splited_str[0];
    genome    = splited_str[1];
    
    //Matrix<MICell, 2> matrix(10000, 10000);
    MICell* matrix = new MICell[(10000+1)*(10000+1)];

    while(cin) {
        getline(cin, input_line);
        if (input_line.size() == 0) {
            return 0;
        }
        splited_str = split(input_line, '\t');
        if (splited_str.size() != 3) {
            cerr << "Bad gene line!" << endl;
            return 1;
        }
        gene_id     = splited_str[0];
        start_pos   = str2int(splited_str[1]) - 1;
        gene        = splited_str[2];
        length      = gene.size();

        if ((gene_id == prev_gene_id) and (prev_score >= MIN_GOOD_SCORE)) {
            // уже получили значимый score для пары этих gene-genome в предыдущих итерациях
            continue;
        }

        if ((length >=100) and (length <= 10000)) {
            genome_subseq = genome.substr(start_pos, length);
            if (use_external_interpolation) {
                // create & adjust matrices
                //cout << "|" << gene << "|\n";
                //cout << "|" << genome_subseq << "|\n";
                CharWeights gene_weights   = create_char_weights(gene);
                CharWeights genome_weights = create_char_weights(genome_subseq);
                //print_weights(gene_weights);
                //print_weights(genome_weights);
                //cout << "\n";
                adjust_weights(gene_weights);
                adjust_weights(genome_weights);
                //print_weights(gene_weights);
                //print_weights(genome_weights);
                // get koeffs
                std::string db_key = create_db_key_pupy(genome_weights, gene_weights);
                if (interpolation_data.count(db_key) == 1) {
                    raw_score = _nw_score_restricted(genome_subseq, gene, matrix);
                    score = interpolation_data[db_key].calculate(length, raw_score);
                }
                else {
                    std::cerr << "Cannot find interpolation db_key: \"" << db_key << "\"\n";
                    return 1;
                }
            }
            else {
                score = nw_score_mk(genome_subseq, gene, matrix, 10);
            }
            cout << genome_id << "\t" << gene_id << "\t" << score << endl;
            prev_gene_id = gene_id;
            prev_score = score;
        };
    }

    return 0;

}
